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Abstract 

While  the  use  of  relative  orbit  determination  has  made  in  reducing  minimizing 
the  difficulties  inherent  in  tracking  geostationary  satellites  that  are  in  close  proximity, 
the  problem  is  often  compounded  by  stationkeeping  operations  or  unexpected  maneu¬ 
vers.  If  a  maneuver  occurs,  observations  will  no  longer  fit  predicted  data,  increasing 
the  risk  of  misidentification  and  cross-tagging. 

The  goal  of  this  research  was  to  develop  a  model  that  will  estimate  the  mag¬ 
nitude,  direction,  and  time  of  a  suspected  maneuver  performed  by  a  collocated  geo¬ 
stationary  satellite.  Relative  motion  was  modelled  using  Hill’s  equations,  and  least 
squares  estimation  was  employed  to  create  both  a  linear  non-maneuver  model  and 
non-linear  maneuver  model.  Two  sets  of  data  for  an  actual  satellite  collocation  were 
obtained  from  the  Air  Force  Maui  Optical  and  Supercomputing  (AMOS)  site,  con¬ 
sisting  of  differential  right  ascension  and  declination.  Studies  conducted  with  these 
observations,  along  with  simulation  studies,  indicate  that  it  is  indeed  possible  to  per¬ 
form  maneuver  estimation.  It  was  found,  however,  that  the  amount  of  data  required 
for  successful  convergence  is  much  greater  than  that  typically  obtained  for  tracking 
purposes. 
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Maneuver  Estimation  Model  for  Relative  Orbit 


Determination 

I.  Introduction 

1.1  Background 

The  concept  of  the  geosynchronous  orbit,  and  it’s  more  specific  counterpart,  the 
geostationary  orbit,  has  been  around  for  more  than  a  century.  While  Arthur  C.  Clark 
became  widely  known  for  this  concept  in  October  of  1945,  it  was  first  proposed  in 
the  early  1900’s  by  Russian  theorist  Konstantin  Tsiolkovsky.  This  concept  became  a 
reality  in  July  of  1963,  when  Syncom  2  became  the  first  operational  geosynchronous 
communications  satellite  [7]. 

Since  then,  the  geosynchronous  orbit  (GEO)  regime  has  proven  to  be  an  invalu¬ 
able  asset,  so  much  so  that  the  number  of  satellites  placed  in  these  orbits  has  risen 
dramatically.  While  demand  for  this  capability  continues  to  increase,  the  available 
number  of  orbital  slot  allocations  continues  to  decrease.  Consequently,  many  organi¬ 
zations  are  choosing  to  collocate  satellites  in  the  same  slot.  In  addition  to  intentional 
collocation,  cases  now  exist  where  satellites  unwittingly  have  been  placed  in  a  position 
where  one  stationkeeping  box  overlaps  another,  leading  to  an  increased  vulnerability 
of  unintentional  close  approaches  [4].  And  of  course,  as  in  any  orbital  regime,  an  in¬ 
crease  in  space  debris  and  malfunctioning  vehicles  lead  to  another  potential  for  close 
approaches  and  collisions. 

Whether  intentional  or  unintentional,  collocation  and  close  approaches  increase 
the  difficulty  of  identifying  individual  satellites  within  clusters  and  create  the  potential 
for  misidentiheation  and  cross-tagging.  While  various  identification  methods  exist, 
increased  orbit  determination  accuracy  is  a  valuable  way  to  monitor  the  extent  of  close 
approaches,  thus  minimizing  the  need  for  expensive  precautionary  collision  avoidance 
maneuvers.  Relative  motion  has  emerged  as  a  potential  asset  in  supplying  spacecraft 
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identification  information.  By  using  relative  metric  data  from  optical  sensors  and 
relative  equations  of  motion,  spacecraft  separations  can  be  estimated  and  predicted  [5]. 

1.2  Problem  Statement 

While  relative  orbit  determination  has  made  improvements  in  minimizing  the 
inherent  difficulties  in  tracking  objects  in  close  proximity,  the  problem  is  often  com¬ 
pounded  by  stationkeeping  operations  or  unexpected  maneuvers  performed  by  one 
or  more  satellites.  If  a  maneuver  takes  place,  the  observations  will  no  longer  fit  the 
predicted  data  and  misidentihcation  and  cross  tagging  are  problems  once  again. 

1.3  Method  of  Solution 

The  goal  of  this  research  is  to  create  a  model  that  will  estimate  the  magnitude, 
direction,  and  time  of  a  suspected  maneuver  performed  by  a  collocated  satellite.  Rel¬ 
ative  orbit  determination  and  least  squares  estimation  are  employed  to  create  both  a 
linear  non-maneuver  model  and  non-linear  maneuver  model.  Observations,  obtained 
from  the  Air  Force  Maui  Optical  and  Supercomputing  (AMOS)  site,  consist  of  dif¬ 
ferential  right  ascension  and  declination,  and  dynamics  will  be  modelled  using  Hill’s 
equations. 
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II.  Literature  Review 


2.1  Optical  Observations 


2.1.1  Fundamentals  -  Celestial  Sphere  Geometry.  [8,15] 

Most  celestial  objects  observed  from  Earth  are  at  a  distance  many  times  greater  than 
that  of  the  earth’s  radius,  giving  each  an  apparent  fixed  position  on  the  inner  surface 
of  a  celestial  sphere.  This  discussion  then,  will  begin  with  definitions  associated  with 
the  celestial  sphere,  see  Figure  2.1. 


Figure  2.1:  Geometry  of  the  Celestial  Sphere 


The  celestial  equator  is  the  plane  passing  through  the  earth’s  equator  which 
intersects  the  celestial  sphere.  The  celestial  poles  are  defined  as  the  intersection  of 
the  celestial  sphere  with  the  rotation  axes  of  the  earth,  both  north  and  south.  A  great 
circle  is  the  intersection  of  the  celestial  sphere  with  any  plane  passing  through  the 
center  of  the  sphere.  An  hour  circle  is  one  such  great  circle.  Hour  circles  are  defined 
as  great  circles  that  pass  through  the  celestial  poles  and  are  thus  perpendicular  to  the 
celestial  equator. 

Another  important  concept  involves  the  revolution  of  the  earth  about  the  sun, 
or  as  seen  from  the  observer,  the  apparent  motion  of  the  sun  about  the  earth.  The 
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mean  plane  of  the  earth’s  orbit  around  the  sun  is  called  the  ecliptic  plane  and  is  in¬ 
clined  approximately  23.5  deg  with  respect  to  the  earth’s  mean  equator.  The  apparent 
motion  of  the  sun  on  the  celestial  sphere  as  seen  from  the  earth  follows  the  ecliptic 
plane.  The  intersection  of  the  earth’s  equatorial  plane  with  the  ecliptic  plane  creates 
a  line  connecting  the  equinoxes,  or  line  of  nodes.  The  point  of  intersection  where 
the  apparent  sun  crosses  the  celestial  equator  from  south  to  north  is  termed  the  first 
point  of  Aries,  or  the  vernal  equinox. 

This  framework  provides  a  method  for  defining  the  position  of  objects  in  space. 
Using  the  celestial  equator  as  the  fundamental  plane  and  the  vernal  equinox  as  a 
reference  point  or  principal  direction,  it  is  possible  to  define  two  angular  coordinates 
which  uniquely  determine  the  direction  of  an  object  with  respect  to  the  celestial 
sphere.  Two  such  angular  coordinates  are  right  ascension  and  declination.  Right 
ascension  is  the  angle  measured  east  from  the  vernal  equinox  to  the  particular  hour 
circle  passing  through  the  object  being  observed.  Declination  is  the  angle  measured 
from  the  celestial  equator  to  the  position  of  the  object.  [8, 15] 

2.1.2  The  Raven  Telescope.  The  Raven  telescope,  developed  by  the  Air 
Force  in  1995,  is  an  optical  sensor  designed  to  provide  high  accuracy,  deep  space 
observations.  Made  entirely  of  commercial-off-the-shelf  (COTS)  products  and  fully 
automated,  the  .36m  Raven  performs  ballistic  tracking  with  subarcsecond  accuracy. 
The  telescope  captures  images  using  a  charge-coupled  device  (CCD)  with  a  field  of 
view  of  43  x  29  arcseconds.  Given  this  field  of  view,  it  is  possible  to  simultaneously 
track  multiple  satellites  in  geosynchronous  orbits  [12].  CCD  images  of  satellite  clusters 
provide  more  accurate  metrics  of  vehicle  separation  since  error  sources  introduced  in 
the  observation  process  are  common  to  each  satellite  in  that  frame.  Satellite  position 
is  then  computed  using  astrometry  techniques. 

Originally,  only  one  track  mode,  called  sidereal  mode,  was  used.  By  slewing  the 
telescope  at  the  sidereal  rate,  stars  would  appear  as  points  while  satellites  appeared  as 
streaks.  The  end  points  of  each  streak  were  then  compared  with  the  stellar  background 
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and  the  time  of  camera  shutter  opening  or  closing  was  recorded.  While  this  method 
allowed  for  accurate  fits  to  the  star  field,  satellite  endpoint  streak  detection  introduced 
a  lot  of  uncertainty.  The  next  mode  employed  was  called  stare  mode  and  involved 
maintaining  a  stationary  telescope  position  for  the  duration  of  the  image.  In  this 
mode,  all  of  the  stars  as  well  as  the  non-geostationary  satellites  appeared  as  streaks, 
resulting  in  undesired  image  clutter.  In  addition  to  this  issue,  stare  mode  did  not 
alleviate  the  streak  endpoint  detection  uncertainty.  The  most  recent  tracking  mode 
development  successfully  employed  by  Raven  is  termed  ballistic  or  rate  track  mode.  In 
this  mode,  the  telescope  follows  the  satellite  for  the  duration  of  the  image,  producing 
a  point  for  the  satellite  and  streaks  for  stars,  thus  replacing  the  endpoint  detection 
issue  with  a  centroiding  approach  [3,12],  A  Raven  image  obtained  using  this  method 
is  shown  in  Figure  2.2. 


Figure  2.2:  Raven  Image 


2.1.3  Topocentric  to  Geocentric  Conversion.  When  using  optical  observa¬ 
tions  it  is  essential  to  note  the  coordinate  system  in  which  they  are  expressed.  A 
geocentric  coordinate  system  has  its  origin  located  at  the  center  of  the  earth  while  a 
topocentric  coordinate  system  has  its  origin  translated  from  the  center  of  the  earth 
to  the  position  of  the  telescope  located  on  the  surface  of  the  earth.  In  addition  to  the 
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translation  of  origins,  the  transformation  between  the  two  frames  involves  rotations 
through  the  local  sidereal  time  and  the  compliment  of  the  geodetic  latitude.  Raven, 
along  with  many  other  earth  based  optical  telescopes,  outputs  observations  expressed 
in  the  topocentric  frame.  Geocentric  observations  are  desirable  for  the  purposes  of  this 
research;  consequently,  it  is  necessary  to  perform  the  appropriate  coordinate  trans¬ 
formation.  The  following  formulation,  shown  in  Vallado  (taken  from  Orbital  Motion 
by  Archie  E.  Roy),  determines  the  geocentric  values  of  right  ascension  (a)  and  dec¬ 
lination  (d)  from  the  same  values  in  the  topocentric  frame,  at  and  dt,  respectively. 
This  formulation  also  requires  the  site  position  magnitude,  rSite, and  the  slant  range, 
p.  For  the  purposes  of  this  research,  p  will  be  defined  using  an  average  range  for  a 
geostationary  satellite  of  39149  km. 

77  cos  (j)gc  sin  at  -  6LST  .  N 

tan  at  —  a= - -1 - — - - - - -  (2.1) 

cos  dt  +  cos  <pgc  cos  at  -  0LST 


where  (j>gc  is  the  geocentric  latitude,  measured  positively  north  from  the  equator,  and 
Olst  is  the  local  mean  sidereal  time,  measured  positively  to  the  east  from  the  site. 
The  temporary  variable  7  is  used  to  complete  the  formulation: 


tan  7 


tan  <f>gc  cos  (^2^) 
cos  (7a  -  e^r) 


(2.2) 


r^iiSL  sin  6ac  sin  dt  —  7 

tan  dt  —  d  =  — - - - 7 — - - - 

gm  cy  7.  _sMe  gm  CQS  ^  _  7 


(2.3) 


For  a  more  rigorous  derivation  of  the  above  equations,  see  Roy(1988, 64-67)  [15]. 


2.2  Least  Squares  Estimation 

The  motion  of  an  orbiting  body  about  the  Earth  is  modeled  using  equations  of 
motion,  the  most  basic  of  which  is  the  two-body  equation.  This  equation  describes 
an  orbit  using  six  orbital  elements.  Hence,  when  determining  the  path  of  an  orbiting 
body  using  optical  measurements  of  right  ascension  and  declination,  it  is  necessary 
to  obtain  at  least  three  measurements  -  six  known  values  result  in  six  solvable  values. 
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Typically  however,  more  than  three  observations  are  obtained,  and  all  are  assumed 
to  contain  some  error.  This  first  became  an  issue  in  the  early  1800’s.  In  response  to 
a  mysterious  data  arc,  later  shown  by  Gauss’  methods  to  be  the  astroid  Ceres,  Gauss 
developed  the  theory  of  probability,  leading  to  the  Principle  of  Maximum  Likelihood 
and  the  method  of  Least  Squares  in  its  full,  non-linear  form  [17]. 


2.2.1  The  Principle  of  Maximum  Likelihood.  [17] 

Given  N  independent  measurements  ay  of  the  value  we  want  to  know  Xo,  the 
joint  probability  is  simply  the  product  of  the  individual  Gaussian  distributions 


P(X!,X  2,  ...XN) 
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where  cq  is  the  standard  deviation  of  each  instrument  used  to  obtain  measurements. 
Gauss  then  proposed  that  since  the  true  value  xo  is  unobtainable,  it  is  reasonable  to 
replace  xq  with  an  estimate  x. 
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Subsequently,  the  closest  value  to  the  truth  is  found  when  the  estimate  maximizes  the 
probability  of  having  gotten  the  true  value.  This  is  shown  mathematically  when  the 
argument  of  the  exponential  in  Equation(2.5)  is  minimized,  or  its  derivative  is  driven 
to  zero. 
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Thus  the  name  the  method  of  least  squares. 


(2.6) 


2.2.2  Linearized  Dynamics  and  the  State  Transition  Matrix.  [17] 

In  most  estimation  problems,  the  estimate  of  multiple  values  is  required.  One  portion 
of  this  thesis,  for  example,  is  interested  in  determining  the  position  and  velocity 
differences  between  two  satellites.  These  values  are  typically  organized  in  the  form  of 
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a  system  state  vector,  x.  How  the  state  vector  changes  with  time  can  be  expressed 
using  the  equations  of  motion 

dx  . 

*=»<*•*>  (2'7) 

or  the  explicit  solution  in  terms  of  the  initial  state  and  time 

x(t)  =  h(x(t0),t)  (2.8) 


Either  of  these  equations  specify  how  the  true  state  x0,  the  estimated  state  x,  and 
any  nearby  trajectories  change  with  time. 

Assuming  that  the  estimate  of  the  true  state  is  close  to  the  actual  state,  it  is 
helpful  to  determine  how  two  close  orbits  behave  with  respect  to  each  other.  This 
allows  x  =  x0  +  Sx  to  be  substituted  into  the  equations  of  motion.  Equation  (2.7) 
then  becomes 

—  =  g(x0  +  Sx,t)  (2.9) 

Expanding  g  in  a  Taylor  series  about  the  true  trajectory  yields 

fix 

—  PS  g(x0,  t)  +  V2:g(x0,  i)5x  +  0(2)  (2.10) 

leading  to  the  well  known  Equations  of  Variation,  expressed  as 

— <5x  =  A(t)5x.  (2-11) 

where 

A(t)  =  VxgL(t)  (2.12) 

Because  the  Equations  of  Variation  are  linear  ordinary  differential  equations,  a 
solution  can  be  expressed  as  the  sum  of  the  components  Sxi  of  hx  at  to  multiplied  by 
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each  individual  solution  vector  function  fa: 

N 

Sx(t)  =  y^5xi(to)fa(t)  (2.13) 

i= 1 

Equation  (2.13)  can  now  be  simplified  by  combining  the  individual  components  Sxi 
into  the  vector  hx(t0)  and  the  associated  individual  solutions  <pi  into  the  matrix  $ 
yielding 

dx(t)  =  ®(t,t0)6x(t0)  (2.14) 

This  equation  then  shows  that  the  state  transition  matrix  $  describes  how  a  change 
in  the  initial  conditions  of  the  state  propagates  forward  to  the  value  of  the  state  at 
some  time  in  the  future.  As  stated  above,  the  system  dynamics  can  also  be  given  in 
the  form  of  an  actual  solution  in  terms  of  the  initial  state  and  the  time,  Equation 
(2.8).  In  this  case,  $  can  also  be  expressed  as 

=  Vx(to)h(x(t0),t)  (2.15) 

2.2.3  Linear  Least  Squares.  [17] 

When  estimating  the  linear  state  of  a  system  at  the  epoch  time  to  it  is  necessary  to 
first  look  at  the  observations  z *(£$)  taken  at  each  observation  time  tj.  It  is  assumed 
that  each  observation  vector  z,  is  independent  of  all  other  observation  vectors  and  has 
an  associated  instrumental  covariance  Qi  measuring  the  degree  of  this  independence. 
Assuming  also  that  there  is  a  linear  relationship  between  the  system  state  and  the 
observations  at  any  time  tl:  an  observation  could  be  expressed  with  the  following 
equation,  called  the  observation  relation 

z  i(ti)  =  Hjx(tj)  +  e*  (2.16) 
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where  e*  is  the  true  error  in  the  observations.  It  is  then  possible  to  insert  the  system 
dynamics  into  the  observation  relation 


Zi(fi)  =  #>*(*>,  fo)x(fo)  +  e. 


(2.17) 


Simplifying  further, 


z  i(ti)  =  Tjx(to)  +  ej 


(2.18) 


where  T;  =  #*$(£;,  t0)- 

It  is  then  common  to  assemble  all  vectors  and  matrices  for  the  N  measurements 
into  larger  matrices. 
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(2.19) 


Using  the  method  followed  in  Wiesel  [17],  the  estimate  of  the  state  vector  at  the  epoch 
time  x(£o)  can  be  stated  as  follows 


x(t0)  =  (T1  Q-'Ty'T1  Q-1  z 


(2.20) 
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with  the  measure  of  accuracy  coming  from  the  covariance  matrix  expressed  as 


Pi(t  „)  =  (T'tQ-'T)-1 


(2.21) 


2.2.4  Nonlinear  Least  Squares.  [17] 

Most  problems  in  the  real  world  are  in  actuality  nonlinear  problems,  either  in  their 
dynamics,  their  observation  relations,  or  both.  Linear  systems  for  these  problems  are 
typically  developed  by  making  certain  assumptions.  While  the  linear  system  of  equa¬ 
tions  may  be  sufficient  for  some  applications,  a  more  rigorous  analysis  is  often  required 
prompting  the  usage  of  nonlinear  dynamics  and  nonlinear  observation  geometry. 

An  explicit  solution  for  the  system  dynamics  can  be  represented  as 

x(t)  =  /i(x(f0),f)  (2.22) 

As  Wiescl  explains,  assuming  the  dynamics  are  deterministic,  it  should  follow  that 
their  linearization  about  a  reference  trajectory  xrey 

<5x(t)  =  $(£,  to)Sx(to)  (2.23) 

is  valid,  where  hx  is  the  desired  change  in  the  reference  trajectory  that  will  make 
the  reference  trajectory  equal  to  the  true  trajectory.  Since  the  true  trajectory  x0  is 
unobtainable,  5x  actually  corrects  the  reference  trajectory  into  the  closest  possible 
estimate  of  the  true  trajectory.  Also,  the  state  transition  matrix  <f>,  is  defined  as  the 
gradient  of  the  solution  with  respect  to  the  initial  conditions 


<!>(£,  t0)  =  VXtQ  h(x(t  0) ,  t)  (2.24) 

The  observation  relation,  a  function  that  will  predict  the  observations  given  the  state 
vector,  can  be  expressed  as 

Zi(ti)  =  G(x(ti),ti)  (2.25) 
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where  z,(ij)  are  the  measurements  taken  at  different  observation  times  t,t. 

In  any  measurement  a  certain  amount  of  error  will  be  present.  The  value  z0  is 
the  vector  containing  perfect  measurements  which  would  give  the  true  state  Xo  and  z 
is  the  vector  of  actual  measurements  which  would  give  the  imperfect  observed  state 
x.  Assuming  that  the  true  error  in  the  data  goes  to  zero  as  the  true  error  in  the  state 
goes  to  zero,  the  true  error  in  the  actual  data  can  be  represented  as 

e  =  z  —  z0  =  G(x,  t)  —  G(x0,  t) 

=  G(x0  +  Sx,  t)  -  G(x0,t)  (2.26) 

dGc  .  . 

where  x  =  Xo  +  hx  and  the  last  line  in  this  set  of  equations  relates  the  error  in  the 
state  to  the  error  in  the  reference  trajectory.  Since  it  is  assumed  that  the  residual  r 
will  approximate  the  true  error  e,  the  equation  for  the  residual  becomes 

Yi  =  z i  -  G (Xref(ti),ti)  (2.27) 

H  is  then  defined  as 

dG 

Hi  =  —(xref(ti),ti)  (2.28) 

Using  the  same  form  as  the  last  line  in  Equation  (2.26),  recognizing  that  the  residual 
is  linearly  related  to  Sx,  and  recalling  Equation  (2.23),  the  equation  for  the  residual 
then  becomes 


r*  «  Hi5x(ti)  =  Hi$(ti,t0)6x(t  0) 


(2.29) 


=  Ti5x(t0) 


The  final  results,  being  in  the  same  form  as  the  linear  least  squares  case,  can 
then  be  written  as 

5x(t0)  =  (TTg-1T)"1TTQ”1r  (2.30) 
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(2.31) 


Pi*  =  (TtQ-'T)-1 

and  the  final  estimate  of  the  trajectory  is  then 

x(t0)  =  xre/(f0)  +  MO) 


(2.32) 


2.3  Relative  Motion 

Relative  motion  describes  the  position  of  one  satellite  with  respect  to  another. 
With  more  satellites  placed  in  the  GEO  belt,  relative  motion  has  become  an  increas¬ 
ingly  important  tool  for  analysis  of  orbit  determination  and  satellite  position. 

2.3.1  Equations  of  Motion.  Relative  equations  of  motion  were  developed 
in  1878  when  Hill  derived  a  set  of  equations  to  describe  the  moon’s  orbit  around 
the  earth.  These  equations  were  then  modified  by  Clohessy  and  Wiltshire  in  1960 
to  describe  relative  motion  in  rendezvous  operations.  With  the  aid  of  several  as¬ 
sumptions,  these  equations  are  developed  based  on  position  and  velocity  differences 
between  two  objects  and  can  be  solved  analytically.  The  assumptions  associated  with 
this  formulation  are  as  follows 

1.  The  reference  orbit  is  circular 

2.  Earth  is  spherically  symmetric 

3.  The  distance  between  objects  is  close  when  compared  to  their  orbital  radii 

Due  to  the  limiting  assumptions  associated  with  Hill’s  equations,  many  other 
sets  of  equations  have  been  developed.  One  such  set  is  know  as  the  Cluster  Orbits 
With  Perturbations  Of  Keplerian  Elements  (COWPOKE)  equations.  Using  mean 
Keplerian  elements  and  element  differences,  a  method  of  expressing  the  relative  equa¬ 
tions  of  motion  for  objects  in  non-circular  orbits  has  been  developed  [13]. 

2.3.2  Differential  Orbital  Element  Effects.  [10] 

When  considering  relative  motion  between  two  or  more  satellites,  an  understand- 


2-11 


ing  of  the  effects  of  differential  orbital  elements  is  helpful.  The  first  consideration 
in  maintaining  a  formation  of  satellites  is  the  semimajor  axis  a.  Orbits  with  differ¬ 
ent  semimajor  axes  have  different  periods,  resulting  in  rapid  formation  dispersion. 
Avoiding  this  dispersion  requires  that  the  mean  semimajor  axis  for  each  satellite  be 
the  same.  Next,  inclination  i  differences  result  in  out-of-plane  separation  at  higher 
latitudes.  Because  each  orbit  then  passes  over  different  portions  of  the  earth,  the 
effects  of  J2  on  each  orbit  will  be  slightly  different  as  well,  resulting  in  differences  in 
nodal  precession  rates.  These  differences  cause  orbital  plane  separation  resulting  in 
increased  formation  separation.  Differences  in  right  ascension  of  the  ascending  node 
f!  result  in  maximum  satellite  separation  at  the  equator.  Coplanar  satellites  achieve 
along-track  separation  by  differences  in  mean  anomaly  M,  while  satellites  in  eccentric 
orbits  can  achieve  radial  separation  by  differences  in  the  argument  of  perigee  u. 

2.4  Relative  Motion  Applications  in  GEO 

As  organizations  become  aware  of  the  slot  allocation  problem  in  GEO,  more 
research  is  being  done  on  existing  satellite  clusters  and  unintentional  close  approach 
encounters.  This  section  will  touch  on  a  few  such  research  efforts. 

2-4-1  GEO  Clusters.  As  technology  develops  and  demand  increases,  avail¬ 
able  slots  in  the  geosynchronous  belt  become  more  limited  in  number.  Organizations 
involved  in  communications  who  own  one  of  the  coveted  slots  look  for  ways  to  best 
utilize  and  exploit  the  limited  space  they  own.  Consequently,  more  organizations 
choose  to  collocate  multiple  satellites  along  the  same  longitude.  EUTELSAT  is  one 
such  organization  [9].  To  ensure  consistent  radio  and  television  broadcasting  and 
increase  multi-mission  capabilities,  EUTELSAT  has  concentrated  five  satellites  on 
15  deg  East.  Given  the  system  requirements,  those  associated  with  the  orbit  control 
strategy  provided  the  driving  constraints.  Two  categories  of  constraints  were  identi¬ 
fied.  The  first  constraint,  dealing  with  constraints  due  to  the  quality  of  RF  service 
provided,  was  easily  simplified  by  coordinating  the  frequency  plans  of  the  satellites 
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within  the  cluster.  The  second,  and  driving,  constraint  dealt  with  those  associated 
with  orbit  control. 

2. 4-1.1  Orbit  Control  and  Stationkeeping  Strategies.  Several  orbit  con¬ 
trol  and  formationkeeping  strategies  have  been  presented  [3,9].  In  deciding  on  which 
strategy  to  employ,  each  organization  must  consider  the  constraints  of  their  system. 
As  in  the  EUTELSAT  system  of  vehicles,  these  constraints  may  include  maintaining 
control  box  parameters,  propellant  consumption,  minimizing  station  keeping  maneu¬ 
vers,  avoiding  simultaneous  maneuvers  of  satellites  within  the  cluster,  and  addition 
or  removal  of  satellites.  Many  commercial  communication  satellites  require  a  control 
box  of  ±0.1  deg  in  both  east/west  and  north/south  directions. 

Longitude  Separation:  This  strategy  involves  separating  the  satellites  by  mean 
longitude.  For  some  applications,  this  technique  is  adequate.  Sauer  proposed  this 
to  be  the  simplest  way  of  collocating  two  satellites  with  independent  missions.  By 
partitioning  a  stationkeeping  deadband  of  ±0.1  deg  into  two  smaller  deadbands,  one  to 
two  additional  satellites  can  be  successfully  collocated  with  an  existing  satellite  with 
minimal  impact  on  existing  stationkeeping  operations  [14].  Pattinson,  on  the  other 
hand,  finds  this  technique  inadequately  susceptible  to  close  approaches  for  clusters 
consisting  of  several  satellites  with  a  stationkeeping  cycle  time  of  two  weeks,  as  with 
the  ELITELSAT  cluster  located  at  15  deg  East  [9]. 

Eccentricity  Vector  Separation:  The  eccentricity  vector,  as  shown  in  Vallado, 
can  be  mathematically  described  as, 

(v2  —  -)  r  —  (r  •  v)v 

e  =  - - — - 

h 

This  vector  has  a  magnitude  equal  to  the  orbit  eccentricity  and  always  points  from  the 
center  of  the  Earth  to  the  orbit  perigee  [15].  Eccentricity  vector  separation  sets  the 
eccentricity  vectors  of  collocated  satellites  in  different  directions  while  maintaining  the 
same  longitude.  This  strategy  produces  radial  separation  over  the  orbit  which  may 


(2.33) 
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seem  adequate,  however,  several  problems  may  occur.  When  the  satellites  are  sepa¬ 
rated  in  the  radial  direction  alone  obstruction  of  one  satellite  by  the  other  is  possible, 
leading  to  a  disruption  of  mission  operations.  In  addition  to  this,  significant  position 
errors  may  develop  due  to  orbit  determination  accuracies  and  maneuver  performance 
uncertainties  which  would  then  lead  to  unacceptable  error  ellipsoids. 

Combined  Eccentricity  and  Inclination  Vector  Separation:  To  deal  with  the  prob¬ 
lems  created  by  eccentricity  vector  only  separation,  a  difference  in  inclination  is  com¬ 
monly  added  to  the  satellite  cluster.  As  defined  by  Pattinson  [9],  the  direction  of 
the  inclination  vector  is  determined  by  projecting  the  orbit  pole  onto  the  earth’s 
equatorial  plane.  It  is  possible  to  orient  the  vector  separation  such  that  any  occupa¬ 
tions  are  avoided.  This  configuration,  however,  is  undesirable  since  the  risk  of  close 
approaches  is  still  prevalent  when  radial  separation  ceases  to  exist.  Another  configu¬ 
ration  places  the  inclination  vector  perpendicular  to  the  eccentricity  vector,  allowing 
for  satellite  separation  in  the  north/south  direction  when  radial  separation  does  not 
exist.  While  occupations  are  still  possible,  this  configuration  is  preferred  over  the 
previous  configuration  since  close  approach  constraints  are  mandated  while  the  prob¬ 
ability  of  occupations  is  typically  small.  Given  these  two  configurations,  it  is  then 
possible  to  optimize  the  orientation  of  the  eccentricity  and  inclination  vectors  such 
that  the  close  approach  constraint  is  met  while  minimizing  occupations. 

2-4-2  Close  Approaches.  Unintentional  close  approaches  are  becoming  more 
frequent  as  the  GEO  belt  becomes  more  populated.  Slot  allocation,  as  determined 
by  the  International  Telecommunications  Union,  is  based  primarily  on  separation  of 
operational  frequencies  while  physical  proximity  is  often  overlooked.  This  has  placed 
satellites,  operated  by  different  organizations,  in  slots  with  overlapping  longitude 
stationkeeping  boxes  [4],  In  addition  to  slot  allocation,  there  remains  the  ever  present 
possibility  of  vehicle  failure  resulting  in  uncontrollable  drift.  In  August  of  1997, 
an  uncontrolled  communications  satellite,  Telstar  401,  came  within  12  kilometers  of 
GOES-10,  an  operational  meteorological  satellite  used  by  the  National  Oceanic  and 
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Atmospheric  Administration  (NOAA)  [11],  Since  that  time,  Telstar  401  has  had 
over  100  close  approaches  with  operational  satellites,  at  distances  as  close  as  two 
kilometers  [1]. 

2.5  Chapter  Outline 

With  this  information  established  it  is  now  possible  to  continue  with  the  subse¬ 
quent  chapters. 

Chapter  Three  gives  the  overall  methodology  for  the  non-maneuver  and  ma¬ 
neuver  models.  Derivation  of  the  observation  function  is  explained,  and  the  state 
vector  and  dynamics  for  each  model  are  discussed.  The  algorithm  involved  in  each 
estimation  run  is  also  presented. 

Chapter  Four  discusses  the  simulation  and  real  data  sets  used  in  this  research. 
Simulated  data  was  created  for  one  non-maneuver  scenario  and  two  maneuver  scenar¬ 
ios.  In  addition  to  simulated  data,  two  sets  of  data  from  actual  collocated  satellites 
are  evaluated  and  analyzed. 

Chapter  Five  summarizes  the  results  of  this  study  and  offers  recommendations 
for  future  research. 
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III.  Methodology 

This  chapter  will  discuss  the  process  by  which  the  relative  separation  between  two 
satellites  will  be  estimated  along  with  a  possible  maneuver  and  the  maneuver  time. 
A  discussion  on  coordinate  transformations  from  a  body-fixed  frame  to  the  geocentric 
inertial  frame  is  followed  by  the  formulation  of  the  observation  matrix.  Next  the 
system  dynamics  for  both  the  linear  and  non-linear  models  is  discussed.  Finally  the 
algorithm  implemented  in  the  computer  model  is  stated. 


3. 1  Coordinate  Transformation 

Geocentric  inertial  right  ascension  and  declination  are  expressed  in  the  earth 
centered  inertial  frame,  with  unit  vector  nx  pointing  in  the  direction  of  the  vernal 
equinox,  unit  vector  hz  going  through  the  North  Pole,  and  unit  vector  ny  completing 
the  right  hand  rule.  Hill’s  equations,  on  the  other  hand,  are  formulated  using  a 
body  fixed  frame,  with  unit  vector  er  pointing  radially  away  from  the  center  of  the 
earth,  unit  vector  eo  in  the  along  track  direction,  and  unit  vector  ez  in  the  cross  track 
direction. 

In  order  to  properly  formulate  the  observation  function,  consisting  of  differential 
right  ascension  and  differential  declination,  the  position  vectors  of  each  satellite  must 
be  expressed  in  the  same  frame.  The  transformation  between  these  two  frames  consist 
of  rotations  involving  the  right  ascension  of  the  ascending  node  fl,  inclination  i,  and 
the  argument  of  latitude  u,  defined  as  the  angle  measured  between  the  ascending  node 
and  the  satellite  position  vector.  The  rotation  matrix  C  from  the  body-fixed  frame 
to  the  geocentric  inertial  frame  can  be  written  as 
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Performing  the  matrix  operations  gives  the  following  expression  for  C 


C 
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—  sin  u  cos  hi  —  cos  u  cos  i  sin  hi 
sin  i  sin  hi 


cos  u  sin  12  +  sin  u  cos  i  cos  hi 
—  sin  u  sin  hi  +  cos  u  cos  i  cos  12 
—  sin  i  cos  hi 


sin  u  sin  i 

cos  u  sin  i 

cos  i  j 
(3.2) 


3.2  The  Observation  Function  and  Its  Linearization 

in  determining  the  relative  separation  of  two  satellites,  it  is  necessary  to  express 
the  given  optical  data,  right  ascension  a  and  declination  d,  as  a  function  of  the  system 
state,  current  time,  and  observation  geometry.  This  is  achieved  using  the  observation 
function  G.  The  following  discussion  shows  the  formulation  of  G. 

The  position  vector  for  satellite  one  R1;  originating  from  the  center  of  the  earth 
and  expressed  in  the  geocentric  inertial  frame  is 


Ri  —  R\ 


cos  aicosdi  sinaicosdi  sindi 


In  \ 


ri„ 

V*  / 


(3.3) 


where  a  is  the  geocentric  right  ascension  and  d  is  the  geocentric  declination.  See 
Figure  (3.1).  The  position  vector  for  satellite  one  can  also  be  expressed  in  the  body 


Figure  3.1:  Geocentric  Inertial  Frame 
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Sat!  \& 

Figure  3.2:  Body  Fixed  Frame 

frame  as  shown  in  Figure  (3.2).  Before  this  vector  can  be  used,  however,  it  must  first 
be  rotated  to  the  geocentric  inertial  frame 


R-satl  =  (  r0  0  0 


/e,\ 
eg 
\ez  J 


/n,\ 


0  0  )C 


y 

\n*  / 


(3.4) 


Setting  these  two  expressions  equal  to  each  other  yields 


Ri  —  R 


sat  1 


Ri 


^  cos  aq  cos  di  ^ 
sin  an  cos  ah 
y  sin  di  j 


(  r0Cu  ^ 
^oCl2 
y  roh'l  3  J 


(3.5) 


Similarly,  the  expressions  for  satellite  two  are 


R2  —  R2  (  cos  o.2cosd2  sin  a2Cosd2  sin  d2 


y 

n,  / 


(3.6) 


Rsat2  =  (  r0  +  hr  roh#  j 


/  e;  \ 

eg 

\dz  J 


(n\ 


=  (  ro  +  Sr  r0se  6z^j  C 


n,, 


y 

\nz  J 


(3.7) 
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(3.8) 


R2  —  R-sat2 


^  cos  a2  cos  d2  ^ 

(  (r0 

R-2 

sin  a2  cos  d2 

= 

(ro 

y  sin  d2  j 

\  ( ro 

+  6r)Cu  +  vqSOCzi  +  SzC^i 
+  Sr)C\2  +  ro59C22  +  dzC^ 

+  5r)Ci$  +  To5dC23  +  SzC^s  j 


The  differential  position  vector,  and  subsequently  the  differential  right  ascension  and 
declination,  can  then  be  found  by  subtracting  equations  (3.5)  and  (3.8) 


Ri  —  R2 


R 


sat  1 


R 


sat2 


(3.9) 


producing  the  following  three  equations 


_R2(cos  a2  cos  rf2)  —  Ri(cos  aq  cos  d\)  =  (r0  +  Sr)Ci2  +  r056C2i  +  SzC3i  —  R\C\  i  (3.10) 

i?2(sino;2cos(i2)  —  i?i(sina!  cosdi)  =  (r0  +  5r)C12  +  r056C22  +  SzC32  —  i?iC'i2  (3.11) 
i?2  sin  d2  —  R\  sin  d\  =  (tq  +  Sr )C  13  +  Tq86C2^  +  5zC 33  —  R\C\^  (3.12) 

Beginning  with  Equation  (3.12),  it  is  possible  to  obtain  an  expression  for  the 
differential  declination,  8d.  Recalling  the  following  trigonometric  formula  for  some 
values  A  and  B 

sin  A  +  B  =  sin  A  cos  B  +  cos  A  sin  B  (3.13) 

and  setting  R2  =  R\  +  8R  and  d2  =  d±  +  5d,  the  left  side  of  Equation  (3.12)  becomes 

(Ri  +  5R)(sind\  cos  8d  +  sirufdcosdi)  —  R\  sin  d\  (3-14) 
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Substituting  this  in  for  the  left  side  of  Equation  (3.12)  and  using  the  small  angle 
assumption,  the  following  equation  results 


Ri  sin  Gh+sin  d\5R+R\  cosdi5d—R\  sindi  =  RiCis+5RCi3+Ri59C23+5zCs3  —  R1C13 

(3.15) 

where  Ri  =  r0  and  5R  =  5r.  Simplifying  this  equation  yields 


r  r,  ,  5R  5R  5z 

cos  d\dd  =  —  sm  di  — — K  —613  +  06  623  +  —  633 

III  -Ttl  III 


(3.16) 


But,  recalling  Equation  (3.5),  sin d\  =  C13  so  the  5R  terms  cancel  out.  Consequently, 
the  following  expression  for  5d  results 


r  ,  cosusm«„  cosz 
5d  = - : — 59  H - r^-8z 


COS  d1 


cos  d\R\ 


(3.17) 


Applying  the  same  method,  an  expression  for  differential  right  ascension,  5a , 
can  then  be  solved  for  using  Equations  (3.10)  and  (3.11). 


cos  cos  d\5R  —  R\  sin  a\  cos  d\5a  —  R\  cos  a\  sin  d\5d  = 
RiCn  +  8RC\  1  +  R,x59C2\  +  SzC31  -  RxCu 


(3.18) 


sin  cos  d\5R  +  R\  cos  an  cos  d\5a  —  R\  sin  a\  sin  d\5d  = 

RiC12  +  5RC\2  +  R\59C22  +  5zC32  ~  R1C12  (3.19) 

Multiplying  Equation  (3.18)  by  sinai,  Equation  (3.19)  by  cosai,  and  subtracting 
Equation  (3.18)  from  (3.19)  yields 

Ri  cos  d±5ai  = 

( C12  cos  ai  —  C\\  sin  a\)5R  +  ( C22  cos  a\  —  C2\  sin  a\ )R\59  + 

(C32  cos  a\  —  C31  sin  a\)5z  (3.20) 
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Recalling  again  Equation  (3.5), 


cos  aq  cos  d\  =  C\  \ 
sin  aq  cos  d\  =  C\2 

and  noting  these  expressions  in  the  5R  term  of  Equation  (3.20) 


C 12  cos  aq  —  Cn  sin  aq  =  0 


(3.21) 


The  expression  for  differential  right  ascension  then  becomes 


5a  =  — — —  (  (C2 2  cosaq  —  C2i  sinaq)50  +  (C32  cosaq  —  C31  sinaq)—  )  (3.22) 

cos  di  V  R\) 


This  expression  can  be  further  reduced  by  recognizing  again  from  Equation  (3.5) 


C 


cos  ctq 
sin 


11 


cos  d\ 
C12 
cos  d\ 


Considering  the  59  component  of  equation  (3.20) 

n  n  .  C22C11  —  C21C12  cos  i 

C22  cos  ai  -  C2i  sm  oq  =  - - - = - - 

cos  ai  cos  d\ 


(3.23) 


Similarly,  considering  the  z  component  of  Equation  (3.20) 


C32  cos  —  C31  sin  = 


C32C11  -  <^31^12  _  sin  i  cos  u 

cos  diRi 


cos  di 


(3.24) 


Finally,  the  expression  for  differential  right  ascension  becomes 


5a  = 


cos2  di 


...  r/1  sinicosM. 

cos  ( i)59 - - - 5z 

R 1 


(3.25) 
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Having  solved  for  dd  and  da,  the  observation  relation  G  is  now  written  as 


G  = 


^(cos  (i)se-^Sz) 

cos_usinz  'T/3  I  cos i  r 
cos  d\  cos  d\ R\  ^ 


(3.26) 


Due  to  the  linearization  about  satellite  one  using  Taylor  series  expansion  in  the  de¬ 
velopment  of  G,  its  linearization  H  is  simply 


H 


dG 

dX 


0 

0 


cos  i 

R\  cos2  d\ 

cos  u  sin  i 
Ri  cos  di 


sin  i  cos  u 
R\  cos2  d\ 

cos  i 

Ri  cos  d\ 


0  0  0 
0  0  0 


(3.27) 


where  the  state  vector  is  defined  as 


X 


dr  r069  5z 


5r  r069 


(3.28) 


3.3  Linear  System  Dynamics 

The  relative  motion  between  two  co-located  satellites  in  geostationary  orbit  can 
be  described  using  Hill’s  equations.  Assuming  a  circular  reference  orbit  and  a  small 
distance  between  satellites,  Hill’s  equations  take  the  following  form: 

dr  —  2  nr0dd  —  3  n2dr  =  0 

r0d9  +  2  ndr  =  0  (3.29) 

dz  +  n2dz  =  0 
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As  Wiesel  [16]  explains,  these  equations  can  be  solved  analytically 


hr(t) 

59(t) 

5z(t ) 
5r(t) 

59(t) 

5z(t) 


—  I  —  r0590  +  3hr0  ]  cos  nt  +  sin  nt  +  4 hr0  H — r0590 
n  /  n  n 


590  -  (^3S90  +  t  + 

r  hi0  . 

ozo  cos  nt  H - sin  nt 

n 


4S90  6 hr0 


n 


r0 


2 hr0  2  r. 

sin  nt  H - cos  nt - hr0 


nr0 


2r059o  +  3nhr0  )  sin  nt  +  hr0  cos  nt 
f  6  n5r0 


-3  +  l 

r0  )  V  ro 

-hzon  sin  nt  +  Szq  cos  nt 


+  45 9o  )  cos  nt - —  sin  nt 

r0 


nr  o 
(3.30) 


The  solution  can  then  be  put  into  matrix  form,  represented  by  the  6x6  matrix  <&hm- 

*"■  'j  (3.31) 

The  solution  for  each  position  component  (hr,  r059,  5z)  is  found  by  propagating  the 
initial  position  components  (hr0,  roh^o,  h20)  forward  using  the  3x3  matrix  <hrr  and 
the  initial  velocity  components  (hr0,  r0S90,  5zq)  forward  using  <!>„,.  Similarly,  the 
solution  for  each  velocity  component  (hr,  rohd,  Sz)  is  found  by  propagating  the  initial 
position  components  forward  using  and  the  initial  velocity  components  forward 
using  <!>„„.  This  can  be  shown  in  the  following  equations 

hr  =  $rrhr(t  =  0)  +  $rv5v(t  =  0)  (3.32) 

hv  =  $vr5r(t  =  0)  +  $vv5v(t  =  0)  (3.33) 
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Referring  back  to  Equation  (3.31),  the  full  state  transition  matrix  is  written  as 


/ 


4  —  3  cos  0 
6(sin  0  —  0) 


®Hill  = 


0 

3n  sin  0 
6n(cos0  —  1) 


V 


o 


0  0 

1  0 

0  cos  0 

0  0 

0  0 

0  — nsin0 


-  sin  0 

n 

(cos  0  —  1) 
0 

COS  0 

— 2  sin  0 
0 


1(1  -  COS  0) 
-  sin  0  —  40 

o 

2  sin  0 
—3  +  4cos0 
0 


0  ^ 

0 

4  sin  0 

o 

o 

COS  0  j 

(3.34) 


where  n  is  the  mean  motion  of  the  reference  satellite  and  0  =  nt.  The  state  vector 
for  the  linear  model  is  defined  as 


X 


Sr  r0S9  Sz  Sr  r0S6  Sz 


(3.35) 


3.4  The  Maneuver  Model:  Non-linear  System  Dynamics 

If  X0  is  the  initial  state  of  the  system,  as  defined  in  Equation  (3.35),  the  state 
of  the  system  at  any  time  t  prior  to  a  maneuver  Xpre  can  be  determined  using  the 
state  transition  matrix  given  in  Equation  (3.34) 

Xpre  =  R))X0  (3.36) 

The  state  immediately  after  a  maneuver  Npost  can  then  be  defined  as  the  initial  state 
propagated  forward  to  the  time  of  the  maneuver  tm  plus  the  change  in  state  caused 
by  the  maneuver  AX 

post  —  $ Hillitmi  ^o)Xo  +  AX  (3.37) 

Note  that  because  the  satellite  will  not  be  instantaneously  displaced,  AX  will  only 
consist  of  velocity  components.  Using  this  equation,  the  state  of  the  system  X  at 
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some  time  t  after  the  maneuver  can  be  written  as 


X  —  $ Hill{t,tm )  Xpost 

=  $Hill(t,  ^o)Xq  +  tm)AX 


(3.38) 


Simplifying,  X  becomes 


X  —  ^o)Xo  +  &Hui(t,  tm)AX.  (3.39) 

This  can  be  expressed  graphically  as  shown  in  Figure  3.3 


Figure  3.3:  Maneuver  Dynamics 

Because  the  goal  of  this  research  is  to  estimate  the  size  and  direction  of  the 
maneuver  as  well  as  the  maneuver  time,  the  state  vector  in  the  non-linear  maneuver 
model  will  include  the  AX  components  as  well  as  the  change  in  maneuver  time. 

Xm  =  (<5r,  roS0,  8z,  5r ,  ro60,  8z ,  Ary,  Avq,  Avz,  8tm)T  (3.40) 

where  Avr,  Avg,  and  Avz  are  the  changes  in  velocity  in  the  dr,  r^Sd,  and  8z  directions 
respectively,  and  Stm  is  the  change  in  the  estimate  of  the  maneuver  time. 

The  state  transition  matrix  <Fmaneu„er.  then  becomes  a  10  x  10  matrix  with  the 
upper  6x6  being  the  solution  to  Hill’s  equations,  <&huu  as  seen  hi  Equation  (3.34). 
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The  next  6x3  portion  of  $  maneuver  propagates  AX  forward  in  time.  This  will  be 
done  using  the  3x3  velocity  matrices  from  the  solution  to  Hill’s  equations,  <&rv  and 
as  shown  in  Equation  (3.31).  The  final  6  x  1  in  the  upper  portion  of  $  maneuver 
propagates  the  estimate  for  the  maneuver  time.  Recalling  Equation  (3.39)  and  Figure 
3.3,  this  column  is  the  partial  derivative  of  &Hui(t,tm)  with  respect  to  the  maneuver 
time,  which  is  simply  the  partial  derivatives  of  and  <1>OT  (replacing  the  Sr,  rpS9, 
and  Sz  components  with  Avr,  Avg,  and  Avz)  with  respect  to  the  maneuver  time.  The 
equations  are  shown  below 

=  —2r0Avg  cos  [n(t  —  tm)\  —  Avr  cos  [n[t  —  tm)\ 

=  — 4r0Avg  cos  [n(t  —  tm)\  +  2Aiy  sin  [n{t  —  tm)]  +  3r0Avg 
=  —Avz  cos  [n(t  —  tm)}  (3-41) 

=  —2r0Avg  cos  [n{t  —  tm)\  +  A vrn  sin  [n(t  —  tm)\ 

=  4roAu07isin  [n(t  —  tm))  +  2Aiyncos  [n(t  —  tm)] 

=  A vzn  sin  [n(t  —  tm)] 

The  bottom  4x6  (rows  7  through  10,  columns  1  through  6)  are  the  partial  derivatives 
of  the  Av  components  and  the  maneuver  time  with  respect  to  the  position  and  velocity 
components,  which  are  all  zero.  The  final  4x4  (rows  7  through  10,  columns  7  through 
10)  are  the  partial  derivatives  of  the  Av  components  and  the  maneuver  time  with 
respect  to  themselves,  giving  an  identity  matrix.  The  full  &  maneuver  is  shown  below. 


Sr 

Stm 

roS0 

Stm 

Sz 

Stm 

Sr 

Stm 

r0S9 

Stm 

Sz 

Stm 
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$ 


maneuver 


/ 

\ 

<f> 

dr 

A3  rv 

dtm 

$ Hill 

3x3 

$ 

dv 

^  vv 

dtm 

6x6 

3x3 

0 

/ 

^  4x6 

4x4  j 

(3.42) 


3.5  Algorithm 

This  section  will  outline  the  algorithm  that  was  used  in  the  estimation  program 
including  the  type  of  data  used,  initialization  of  the  state  vector  and  the  reference 
orbit,  and  the  least  squares  method. 


3.5.1  Data  Files.  As  the  telescope  tracks  two  satellites,  it  records  a  data  arc 
for  each  vehicle  consisting  of  n  observations  of  right  ascension  and  declination  as  well 
as  the  time  associated  with  each  observation.  Because  relative  equations  of  motion 
solve  for  the  separation  between  two  objects,  it  is  necessary  to  have  the  right  ascension 
and  declination  differences,  5a  and  5d,  respectively,  for  each  time  rather  than  the 
absolute  right  ascension  and  declination  for  each  satellite.  Recalling  Section  2.1.3, 
Raven  observations  are  expressed  in  the  topocentric  frame.  Before  the  observations 
are  differenced,  it  is  necessary  to  convert  the  observations  from  the  topocentric  frame 
to  the  geocentric  frame  as  explained  in  the  same  section.  Once  the  conversion  is 
completed,  the  observations  are  differenced  and  compiled  in  an  n  x  2  matrix  z.  A 
standard  deviation  of  approximately  10  arcseconds,  or  4.848  x  10-5  radians,  is  assigned 
for  each  component  of  the  relative  observations.  These  values  are  then  complied  in 
the  instrumental  covariance  matrix  Q ,  often  referred  to  as  the  observation  weighting 
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matrix 


(4.848  x  1(T5)2 
0 


0 

(4.848  x  10~5)2 


(3.43) 


3.5.2  Initializing  the  State  Vector.  Recalling  Equation  (3.35),  the  state 
vector  for  the  linear  non-maneuver  model  is 


X  =  (Sr,  r0S9, 5z,  Sr, r0S9 , Sz)T 


Because  the  relative  separation  is  assumed  to  be  small,  this  state  is  initialized  using 
a  zero  vector. 

Initialization  for  the  non-linear  maneuver  model  is  less  straight  forward.  Re¬ 
calling  Equation  (3.40),  the  vector  is 


Xm  =  (Sr,  r 059, 5z ,  Sr,  r0S9,  Sz,  Avr,  Ave,  Avz,  5tm)T 


The  first  six  components  are  still  set  with  an  initial  value  of  zero  along  with  Stm.  To 
initialize  the  Av  components,  it  is  necessary  to  recall  the  partial  derivatives  of  position 
and  velocity  with  respect  to  Stm,  as  shown  in  Equation  (3.41).  It  can  be  seen  that 
every  component  in  these  equations  is  dependent  upon  one  of  the  Av  components. 
In  order  to  avoid  observability  problems,  it  then  follows  that  at  least  one  of  the 
three  Av  components  must  be  non-zero.  Referring  to  the  discussion  in  Section  2.4, 
most  stationkeeping  maneuvers  in  GEO  clusters  are  in  the  east-west  or  north-south 
directions,  equating  to  burns  in  the  along-track,  r089,  or  cross-track,  Sz,  directions. 
Various  values  under  2  m/s  are  then  used  as  initial  values  in  either  the  r0S9  or  Sz 
directions.  The  maneuver  time  itself  tm  is  perhaps  the  most  difficult  parameter  to 
initialize.  Without  a  priori  information,  the  entire  estimation  process  must  be  run  for 
different  possible  maneuver  times.  This  is  done  by  looping  through  the  time  vector, 
using  each  observation  time  as  an  initial  value  for  the  maneuver  time. 
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3.5.3  Reference  Orbit.  The  next  step  in  the  algorithm  involves  defining  a 
reference  orbit  for  satellite  one.  This  is  done  using  the  NORAD  Two-Line  Element  set 
(TLE).  Six  parameters  are  extracted  from  the  TLE:  inclination  i,  eccentricity  e,  right 
ascension  of  the  ascending  node  D,  argument  of  perigee  u,  mean  anomaly  M,  and 
mean  motion  n0.  The  semimajor  axis  a  is  solved  for  using  the  mean  motion;  however, 
the  typical  two  body  conversion  cannot  be  used.  This  is  because  the  mean  motion 
given  in  the  TLE  n0  is  actually  the  ’’mean”  mean  motion.  The  following  equations 
recover  a  and  n  from  the  altered  no  given  by  the  TLE: 

=  (~y 

\noJ 

3k2  (3  cos2  io  —  1) 

2ai  (l~e§)3 

=  «i  (i  -  -  %  ~  ffi?)  (3.44) 

3k2  (3  cos2  io  —  1) 

2ao  (l-e^)i 
np 

1  +  ho 

ao 

1  —  ho 

where  ke  =  \/GMe,  with  G  defined  as  the  universal  gravitational  constant  and  Me  as 
the  mass  of  the  earth,  and  k2  =  5.413080  x  10-4  [6]. 

Once  the  complete  set  of  Classical  Orbital  Elements  (COEs)  for  the  reference 
orbit  and  the  mean  motion  have  been  obtained,  it  is  possible  to  continue  with  the  algo¬ 
rithm.  These  values  will  be  used  to  propagate  the  reference  orbit  to  each  observation 
time  and  solve  for  aq  and  d\  of  satellite  one,  as  described  below. 

3.5.4  Least  Squares  Algorithm.  An  observation,  consisting  of  da  and  5d,  and 
its  associated  time  is  read  into  the  loop.  X0  is  propagated  forward  from  the  initial 


a  1 
hi 

ao 

ho 

n 

a 
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time  to  the  observation  time  using  Hill’s  equations  expressed  in  $hm 


X  — 


(3.45) 


The  reference  orbit  is  then  propagated  forward  to  the  observation  time  using  either 
two-body  dynamics  M  =  M0  +  nt0b  or  the  SGP4  propagator.  The  updated  COEs  are 
converted  to  position  r  and  velocity  v  vectors,  which  are  then  used  to  solve  for  the 
right  ascension  and  declination  of  satellite  one,  aq  and  di,  respectively. 


(3.46) 


(3.47) 


These  values  are  then  used,  along  with  i  and  u  (recalling  that  u  is  defined  as  the  angle 
measured  between  the  ascending  node  and  the  satellite  position  vector,  or  u  =  u  +  u 
where  u  is  the  argument  of  perigee  and  v  is  the  true  anomaly),  to  calculate  the 
observation  matrix  G  and  its  linearization  H  as  defined  in  Equations  (3.26)  and 
(3.27).  Once  G  and  H  have  been  obtained  the  matrix  operations  of  the  least  squares 
method,  as  described  in  Section  2.2,  are  performed.  These  steps  are  listed  briefly 
below: 

1.  Solve  for  T,  recalling  T%  =  Hfo 

2.  Calculate  the  residuals,  r*  =  z*  —  G(x) 

3.  Add  new  terms  to  the  sums  of  the  matrix 


and  the  vector 
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This  process  is  repeated  until  each  observation  has  been  processed.  The  remaining 
steps  are  listed  below. 

1.  Calculate  the  covariance  matrix  P$x 

2.  Calculate  the  state  correction  vector  at  the  epoch  time 

P^Y.PQi'  r< 

i 


3.  Calculate  the  new  estimate  of  the  reference  trajectory 

X0(to)  =  X0(fo)  +  SX{tQ) 

This  entire  process  is  then  repeated  for  a  set  number  of  iterations  or  until  the  sys¬ 
tem  has  converged.  Once  this  is  accomplished,  a  new  value  is  used  to  initialize  the 
maneuver  time  and  the  program  is  repeated. 


3-16 


IV.  Simulations  and  Real  Data 


4-1  Relative  Orbit  Determination  Experiment 

In  July  of  2003,  Raven  obtained  images  of  the  DirecTV  4S  and  AMC-4  spacecraft 
collocation  at  101  deg  West  longitude.  Observations  were  taken  during  the  nights 
of  23  -  24  July  and  29  July  -  1  August.  These  observations  were  then  used  by 
researchers  at  the  Air  Force  Maui  Optical  and  Supercomputing  site  to  determine 
if  the  relative  separation  between  two  satellites  could  be  more  accurately  predicted 
than  the  absolute  position  of  each  vehicle  [5].  The  relative  motion  of  the  satellites 
was  estimated  using  the  COWPOKE  equations,  and  the  resulting  differences  in  right 
ascension  and  declination  were  compared  to  the  differences  based  on  the  available  TLE 
from  27  July,  referred  to  as  TLE2.  It  was  determined  that  the  relative  motion  did 
indeed  fit  the  data  better  than  the  solution  provided  by  the  TLE.  One  such  example 
is  shown  in  Figure  4.1. 


Right  Ascension 


Figure  4.1:  COWPOKE  and  TLE  predictions  for  observations  on  30,31  July 

During  the  process  of  conducting  the  study,  however,  it  was  found  that  there  was 
an  occasional  unexplained  shift  away  from  the  actual  position  in  both  the  COWPOKE 
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and  TLE  predictions.  Looking  further  into  the  situation,  it  was  determined  to  be 
highly  likely  that  stationkeeping  maneuvers  had  taken  place;  AMC-4,  referred  to  as 
Sat  1,  may  have  maneuvered  sometime  between  29  and  30  July,  and  DirecTV  4S, 
referred  to  as  Sat  2,  may  have  maneuvered  sometime  between  31  July  and  1  August. 
Figure  4.2  shows  the  differences  in  right  ascension  and  declination,  in  microradians, 
between  the  two  satellites.  A  solution  was  found  using  observations  from  30  and  31 
July  and  propagated  forward  through  1  August.  As  can  be  seen,  the  COWPOKE 
solution  fits  the  data  from  the  30th  and  31st;  however  the  observations  from  the  1st 
have  shifted  away  from  the  predicted  values. 


Right  Ascension 


Figure  4.2:  Fit  to  30f/l  and  31s*,  Predict  to  1st 


4-2  Simulation  Study 

4-2.1  Non-maneuver  Model  Initial  Simulation.  The  first  simulation  created 
was  intended  to  simply  ensure  the  non-maneuver  model  was  working  properly.  The 
reference  orbit  for  Sat  1  was  taken  from  TLE2  while  the  initial  state  vector  Xo, 
consisting  of  the  position  and  velocity  differences  between  the  satellites,  was  chosen 
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to  have  a  separation  of  1750m  in  the  along-track  or  r089  direction. 


X0  = 


0  1750  0  0  0  0 


(4.1) 


The  simulation  was  then  run  for  one  day  (86400s)  with  a  step  size  of  600s.  At  each  time 
step  the  reference  orbit  was  propagated  forward,  and  the  position  and  velocity  vectors 
of  Sat  1  were  computed.  Using  the  method  described  in  Section  3.5.4  and  Equations 
(3.46)  and  (3.47),  the  right  ascension  aq  and  declination  d\  were  then  calculated  for 
Sat  1.  These  values,  along  with  i  and  u,  were  used  to  calculate  the  differential  right 
ascension  and  declination,  8a  and  8cl,  respectively,  using  the  observation  function 
found  in  Equation  (3.26)  and  shown  again  below 


G  = 


=7*  (<=»  (0«  -  ^fc) 

cos  u  sin  i  r/3  i  cos  i  r  ^ 
cos  dp  cosdiRi 


(4.2) 


The  output  of  the  observation  function  was  then  compiled  in  a  hie  and  read  into  the 
non-maneuver,  linear  least  squares  estimation  model. 


Given  the  true  X0  as  shown  in  Equation  (4.1)  the  estimate  converged  on  the 
correct  value  with  near  zero  residuals.  The  root  mean  square  of  the  residuals  ( rrms ) 
shown  below 


T 

r  r 


T  rm.fi 


(4.3) 


was  approximately  4.396  x  10~24.  Figure  4.3  shows  the  trajectory  that  resulted  from 
the  estimated  state  vector,  propagated  forward  for  one  day,  along  with  the  simulated 
observations. 


4-2.2  Maneuver  Model  Initial  Simulation.  There  are  a  number  of  constraints 
inherent  in  dealing  with  operational  optical  systems.  Most  telescopes  obtain  obser¬ 
vations  only  when  the  sky  is  dark  and  the  satellites  are  illuminated  by  the  sun.  This 
effectively  limits  operations  to  night  time  hours.  Weather  is  another  factor  when 
dealing  with  these  systems  since  optical  sensors  cannot  penetrate  cloud  cover.  Due 
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Figure  4.3:  Estimated  and  Simulated  Observations 

to  these  constraints,  a  single  telescope  is  not  capable  of  tracking  a  satellite  or  clus¬ 
ter  of  satellites  over  a  complete  revolution.  This  leads  to  periods  of  time  (ranging 
from  hours  to  days)  where  observations  are  not  available.  While  satellite  motion  is 
deterministic,  maneuvers  can  significantly  alter  a  satellite’s  orbit,  requiring  frequent 
updates  to  the  orbit  determination  prediction.  The  gaps  between  observations  could 
serve  to  inhibit  this  process  by  limiting  the  amount  of  necessary  information  available 
to  accurately  determine  the  current  orbit  of  the  satellites. 

With  this  in  mind,  two  maneuver  simulations  were  created.  The  first  simula¬ 
tion  used  simulated  data  with  perfect  observations  in  similar  quantities  and  times  to 
that  found  in  the  real  DirecTV  4S  and  AMC-4  collocation  data  hies.  Because  two 
maneuvers  were  suspected  between  29  July  and  1  August,  this  study  will  focus  on 
those  four  days.  The  second  simulation  used  continuous  simulated  data  with  perfect 
observations  around  the  entire  orbit  for  the  duration  of  one  day.  Both  simulations 
were  created  with  the  same  initial  state  vector,  as  shown  here 

X0  =  (  0  1750  0  0  0  0  )  (4.4) 
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Another  important  situation  to  explore  is  the  behavior  of  the  maneuver  model 
if  a  maneuver  did  not  take  place.  Hence,  the  final  simulation  involved  evaluating  the 
maneuver  model  when  the  non-maneuver  simulated  data,  described  in  Section  4.2.1, 
was  used. 


4-2.2. 1  Simulated  Non- Continuous  Data.  In  order  to  simulate  the 
actual  observations  of  the  Direct  TV  4S  and  AMC-4  satellite  collocation,  it  was  nec¬ 
essary  to  determine  the  time  and  length  of  each  data  arc  in  the  actual  observation 
files.  Table  4.1  gives  the  approximate  time  (given  in  hour  and  minute  of  that  day) 
and  length  (shown  in  the  number  of  observations)  of  each  arc  for  the  days  of  interest 
as  well  as  the  time  between  observations  within  the  arc,  At. 

Table  4.1:  Data  Arc  Time  and  Length  for  29  July  -  1  August 


July/Aug 

Hour 

Minute 

Obs 

At  (s) 

29 

7 

35-36 

5 

16 

8 

35-36 

5 

16 

9 

32-33 

5 

16 

10 

33-34 

5 

16 

30 

6 

32-33 

4 

16 

6 

45-46 

5 

16 

31 

11 

37-38 

5 

16 

12 

39-40 

5 

16 

13 

31-32 

5 

16 

1 

8 

32-33 

5 

16 

9 

32-33 

5 

16 

10 

34-35 

5 

16 

11 

30-31 

5 

16 

12 

34-35 

5 

16 

13 

35-36 

5 

16 

Three  days  of  simulated  data  were  then  created  to  imitate  the  data  arcs  observed 
in  the  real  data.  This  was  accomplished  using  the  observation  function  G  to  produce 
each  pair  of  observations,  5a  and  5d,  within  the  arc.  Table  4.2  gives  the  time  and 
length  of  each  simulated  arc  as  well  as  the  time  between  observations  within  the  arc. 
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Table  4.2:  Simulated  Data  Arc  Time  and  Length 


Day 

Hour 

Minute 

Obs 

At  (sec) 

1 

7 

35-36 

5 

16 

8 

35-36 

5 

16 

9 

35-36 

5 

16 

10 

35-36 

5 

16 

2 

6 

35-36 

4 

16 

6 

45-46 

5 

16 

3 

11 

37-38 

5 

16 

12 

37-38 

5 

16 

13 

37-38 

5 

16 

As  in  the  real  data  hies,  the  epoch  of  the  estimation  run,  t  —  0,  is  the  time 
of  the  first  observation.  For  instance,  the  epoch  time  of  this  simulation  was  hour  7, 
minute  35  on  the  first  day.  All  times  referred  to  hereafter  will  be  seconds  from  t  =  0 
unless  otherwise  specified  as  the  hour  and  minute  of  a  particular  day. 

A  change  in  velocity  of  2  m/s  was  added  at  t  =  90000  s  (hour  8,  minute  35  on 
Day  2)  to  the  state  X(f  =  90000)  in  the  along-track  or  r056  direction. 

X0  =  X(f  =  90000)  +  (000020)  (4.5) 

Note  that  there  are  no  observations  within  five  hours  of  the  maneuver. 

In  order  to  determine  the  range  of  values  for  the  initial  maneuver  time  guess 
that  will  result  in  convergence  on  the  correct  maneuver  time,  the  estimation  program 
was  looped  through  for  tm  equal  to  each  observation  time.  Once  a  baseline  range  was 
determined,  the  simulation  was  run  again  using  a  smaller  time  step  for  the  initial  tm 
values.  The  state  vector  for  each  run  was  initialized  using  the  truth  state  vector,  that 
used  to  create  the  simulated  data: 

X0  =  (  0  1750  00000200)  (4.6) 
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Results:  An  initial  guess  for  the  maneuver  time  within  approximately  3.5  hours 
before  90000  and  2.5  hours  after  90000  resulted  in  alternating  convergence  on  the 
correct  value  of  tm  =  90000  and  the  incorrect  value  of  tm  =  87138.  The  beginning 
and  end  times  associated  with  each  span  of  initial  maneuver  times,  the  length  of  each 
span,  and  the  resulting  maneuver  time  converged  upon  are  shown  in  Table  4.3. 

Table  4.3:  Convergence  Times  and  Values  for  Simulated  Data  Arcs 


Begin  (s) 

End  (s) 

Span  (min) 

($) 

74520 

74760 

4 

90000 

76200 

76980 

13 

87138 

77040 

77460 

7 

90000 

77520 

80820 

55 

87138 

81660 

82560 

15 

90000 

82620 

87960 

89 

87138 

88020 

91800 

63 

90000 

91860 

95400 

59 

87138 

95460 

95940 

8 

90000 

96420 

98340 

32 

87138 

In  addition  to  the  maneuver  time,  it  is  beneficial  to  look  at  the  solution  of  the 
state  vector.  Table  4.4  gives  the  truth  state  vector  as  well  as  the  solution  to  the  state 
vector  for  each  of  the  above  tm  convergence  values,  with  the  position  components 
given  in  meters,  the  velocity  components  given  in  meters  per  second,  and  the  time 
components  given  in  seconds.  Observing  the  state  vector  solutions,  it  becomes  obvious 
that  tm  =  90000  produces  the  correct  solution.  However,  without  a  priori  knowledge 
of  the  truth  state  vector,  it  would  be  difficult  if  not  impossible  to  determine  which 
solution  is  correct. 

It  is  also  worth  noting  the  standard  deviation,  la,  of  the  maneuver  time  for  each 
solution.  A  maneuver  time  of  90000  produced  a  la  value  of  26257  s  (approximately 
7  hrs  and  20  min)  while  the  tm  =  87138  solution  produced  a  la  value  of  27430  s 
(approximately  7  hrs  and  35  min).  This  essentially  says  there  a  lot  of  uncertainty  in 
these  solutions. 
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Table  4.4:  State  Vector  Solution  for  Convergence  Times 


X0 

Truth 

tm  =  90000 

tm  =  87138 

5r 

0 

-3.0642e-010 

1220.6 

r0Sd 

1750 

1750 

3722.8 

5z 

0 

-1.0123e-014 

-1.1215e-012 

Sr 

0 

7.5276e-015 

-0.030211 

r0S9 

0 

3.7526e-014 

-0.14942 

Sz 

0 

-7.0333e-019 

-2.788e-015 

< 

0 

-8.7235e-014 

-0.79629 

< 

2 

2 

1.9714 

Avz 

0 

7.2749e-019 

2.9755e-015 

Stm 

0 

-1.2407e-009 

-2.2473e-009 

4 -2. 2. 2  Simulated  Continuous  Data.  The  next  simulation  contained 
continuous  simulated  data  with  perfect  observations  (no  noise).  This  simulation  began 
with  the  same  initial  state  vector  X0  as  used  in  the  above  simulations  but  spanned 
one  entire  day  with  observations  every  600  s  (10  min).  The  orbit  was  propagated  for 
half  of  one  day  (43200  s)  and  a  maneuver  of  2  m/s  was  added  at  the  end  of  this  time 
period  ( tm  =  43200  s)  in  the  along-track  or  r0S9  direction. 

X0  =  X(f  =  43200)  +(  0  0  0  0  2  0  )T  (4.7) 

The  orbit  was  then  propagated  for  another  half  day.  Differential  right  ascension  and 
declination  values  were  generated  using  the  observation  function  G. 

The  truth  state  vector  was  used  to  initialize  Xo  in  the  estimation  process,  and 
the  range  of  values  for  the  initial  maneuver  time  guess  that  would  result  in  conver¬ 
gence  on  the  correct  maneuver  time  was  determined,  once  again,  by  looping  through 
the  algorithm  for  tm  equal  to  each  observation  time.  Once  a  baseline  range  was  de¬ 
termined,  the  simulation  was  run  again  using  a  smaller  time  step  between  the  initial 
tm  values. 

Results:  For  the  continuous  data  simulation,  it  was  found  that  the  estimator 
converged  on  the  correct  value  for  tm  if  its  initial  guess  was  within  ±5  hours  of  the 
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actual  time,  43200s.  These  results  produced  a  la  value  of  approximately  11  min  with 
the  solution  to  the  state  vector  shown  in  Table  4.5. 

Table  4.5:  State  Vector  Solution  for  Convergence  Times 


X0 

Truth 

tm  =  43200 

8r 

0 

-1.996e-012 

r059 

1750 

1750 

8z 

0 

2.5015e-016 

Sr 

0 

1.494e-016 

r  086 

0 

2.1308e-016 

8z 

0 

-1.7876e-021 

Avr 

0 

5.039e-016 

Ave 

2 

2 

<1 

0 

6.8903e-020 

Stm 

0 

6.1894e-013 

4 -2. 2. 3  Simulated  Continuous  Data  -  No  Maneuver.  The  final  simu¬ 
lation  involved  running  the  non-maneuver  data  in  the  maneuver  model.  As  stated  in 
Section  4.2.1,  the  initial  state  vector  was  chosen  to  have  a  separation  of  1750m  in  the 
along-track  or  r086  direction. 

X0  =  (  0  1750  0  0  0  0  )  (4.8) 

The  simulation  was  then  run  for  one  day  (86400s)  with  a  step  size  of  600s.  No 
maneuver  was  included  in  this  simulated  data. 

The  state  vector  was  initialized  using  the  truth  position  and  velocity  values, 
that  used  to  create  the  simulated  data.  A  Av  of  2  m/s  in  the  along-track  or  tq60 
direction  was  used  to  initialize  the  change  in  velocity  values.  The  complete  initial 
state  vector  is  shown  in  Equation  (4.9). 

X0  =  (  0  1750  00000200)  (4.9) 
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As  in  the  previous  maneuver  simulations,  the  range  of  values  for  the  initial  maneuver 
time  guess  was  determined  by  looping  through  the  algorithm  for  tm  equal  to  each 
observation  time. 

Results:  For  the  continuous  non-maneuver  data  simulation,  it  was  found  that  for 
each  initial  maneuver  time  that  produced  a  solution,  the  maneuver  model  converged 
on  the  correct  values  for  position,  velocity,  and  Av  with  residuals  on  the  order  of 
1  x  ICR20.  The  model,  however,  was  not  able  to  produce  a  valid  solution  for  the 
maneuver  time.  This  can  be  shown  by  noting  that  the  state  vector  component  Stm 
was  unable  to  converge.  This  makes  sense,  since  there  was  no  maneuver  to  begin 
with.  The  state  vector  for  one  solution  is  shown  in  Table  4.6. 

Table  4.6:  State  Vector  Solution  for  No  Maneuver  Data  Set  in  Maneuver  Model 


X0 

Truth 

tm  =  43200 

5r 

0 

4.9931e-012 

r0S9 

1750 

1750 

5z 

0 

6.8314e-017 

5r 

0 

-1.623e-015 

r056 

0 

-9.9941e-016 

5z 

0 

4.6235e-021 

Aiy 

0 

4.4057e-016 

A  Vg 

0 

1.7448e-016 

Avz 

0 

-1.214e-019 

3tm 

0 

-95371 

4-3  Raven  Data  -  2003 

Recalling  the  relative  orbit  determination  experiment,  discussed  in  Section  4.1, 
two  stationkeeping  maneuvers  were  suspected.  With  a  suspected  maneuver  for  Sat  1 
sometime  between  29  and  30  July  and  another  for  Sat  2  sometime  between  31  July 
and  1  August,  the  observations  were  grouped  into  two  sets  of  data,  each  consisting  of 
three  days  of  observations.  The  first  set  of  data  contained  the  observations  obtained 
for  29-31  July,  while  the  second  set  of  data  contained  the  observations  obtained  for 
30  July  -  1  August. 
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4-3.1  Data  Set  1:  29  -  31  July.  The  first  set  of  data,  consisting  of  observa¬ 
tions  obtained  on  29  -  31  July,  was  used  to  examine  a  possible  maneuver  performed 
by  Sat  1  between  the  29th  and  the  30th  of  July. 

Running  these  observations  through  the  non-maneuver  model,  Figure  4.4  shows 
that  the  linear  solution  to  Hill’s  equations  do  not  give  an  accurate  fit.  Along  with  the 
inexact  fit  to  the  declination,  note  the  curvature  in  the  first  set  of  arcs  shown  in  the 
right  ascension  plot.  The  state  vector  solution  Xq  is  shown  in  Table  4.7. 


Right  Ascension 


Figure  4.4:  Non-maneuver  Fit  to  29  —  31  July 


Table  4.7:  State  Vector  Solution  for  Non-maneuver  Fit  to  29  —  31  July 


X0 

Solution 

dr 

7459.3 

r066 

-62665 

Sz 

5062.6 

Sr 

-0.33898 

r069 

-1.1306 

Sz 

-0.47661 
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Valuable  insight  is  also  gained  by  looking  at  right  ascension  plotted  versus  dec¬ 
lination,  as  shown  in  Figure  4.5.  The  observations  from  29  July  can  be  seen  in  the 


Figure  4.5:  Non-maneuver  Fit  to  29  —  31  July 


upper  left  corner  of  the  plot,  located  at  approximately  -1600  microradians  on  the  right 
ascension  axis  and  150  microradians  on  the  declination  axis  (-1600,150).  Observations 
from  30  July  are  seen  at  approximately  (-1150,-230)  and  those  from  31  July  are  at 
(1050,40).  This  plot  accentuates  the  inaccuracy  of  the  linear  non-maneuver  model  for 
this  data,  noting  particularly  that  the  solution  does  not  fit  the  observations  from  29 
or  30  July.  Along  with  observing  how  well  the  solution  fits  the  data,  it  is  also  possible 
to  examine  various  aspects  of  the  orbital  motion.  Each  ellipse  in  the  solution  plot 
corresponds  to  the  completion  of  one  revolution  of  the  satellites.  The  shift  in  these 
ellipses  are  due  to  the  secular  drift  in  the  along-track  direction  associated  with  each 
orbit.  Note  that  the  size,  orientation,  and  drift  of  each  ellipse  remains  essentially 
constant. 

With  the  discrepancies  in  the  linear  model  apparent,  the  data  hies  are  then  input 
into  the  maneuver  model,  and  the  least  squares  algorithm  is  run  for  initial  maneuver 
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time  values.  These  values  are  obtained  by  stepping  through  the  time  vector  at  60  s 
increments. 

Similar  to  the  simulated  data  presented  in  Section  4.2.2. 1,  convergence  was 
achieved  upon  more  than  one  solution.  The  first  solution  was  tm  =  79477  which 
would  be  during  hour  5  on  30  July.  This  solution  was  converged  upon  for  initial 
maneuver  time  values  of  approximately  7  hours  and  30  min  before  the  assumed  ma¬ 
neuver  and  approximately  2  hours  and  7  minutes  after  the  assumed  maneuver.  The 
model  alternated  converging  on  the  final  two  solutions,  tm  =  125340  (hour  18  on  30 
July)  and  tm  =  168380  (hour  6  on  31  July).  The  times  of  each  convergence  is  shown 
in  Table  4.8  along  with  the  length  of  each  span,  and  the  solution  to  the  state  vector 
for  each  maneuver  time  is  shown  in  Table  4.9. 

Table  4.8:  Convergence  Times  and  Values  for  29-31  July 


Begin  (s) 

End  (s) 

Span  (hr:min) 

tm  (s) 

52400 

87140 

9:39 

79477 

88480 

93060 

1:16 

125340 

93180 

97260 

1:08 

168380 

115420 

120300 

1:21 

168380 

120980 

122840 

0:31 

168380 

122960 

145240 

6:11 

125340 

152360 

182420 

8:21 

168380 

Table  4.9:  State  Vector  Solution  for  Maneuver  Fit  to  29  —  31  July 


X0 

tm  =  79477 

tm  =  125340 

tm  =  168380 

dr 

5380.3 

5410.6 

5529.2 

ro60 

-65297 

-65292 

-65265 

dz 

5673.7 

3903 

3904.3 

dr 

-0.67432 

-0.67883 

-0.69819 

r0dd 

-0.84561 

-0.84922 

-0.86612 

dz 

0.17538 

0.41519 

0.41498 

S- 

< 

0.30420 

-0.2197 

0.72159 

< 

0.028259 

0.094944 

-0.21518 

Avz 

-0.99232 

1.1656 

-1.1655 

dtm 

-3.3343e-007 

-0.00057244 

0.00024853 
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Figure  4.6  shows  how  the  solution  for  tm  =  79477  fits  the  data  versus  time 
along  with  the  TLE  solution  and  the  non-maneuver  solution.  The  maneuver  is  easily 
identified  in  the  declination  plot  at  t  ~  22  hrs. 

The  right  ascension  versus  declination  plot,  as  shown  in  Figure  4.7,  is  also  eval¬ 
uated.  While  this  plot  confirms  that  the  solution  converged  upon  with  the  maneuver 
model  fits  the  data  much  better  than  the  non-maneuver  model,  it  also  reveals  a  con¬ 
siderable  maneuver,  Avz  ~  —1  m/s,  in  the  cross-track,  Sz,  or  N/S  direction.  The 
magnitude  and  direction  of  the  maneuver,  as  seen  in  the  solution  to  the  state  vector 
in  Table  4.9,  could  be  helpful  in  analyzing  and  identifying  a  particular  vehicle. 

The  solution  for  tm  =  125340  is  shown  in  Figure  4.8.  While  the  plots  for  right 
ascension  and  declination  versus  time  look  reasonable,  the  right  ascension  versus 
declination  plot  reveals  that  this  solution  does  not  fit  the  data  as  well  as  the  solution 
from  tm  =  79477,  particularly  in  the  first  set  of  data  arcs.  This  can  be  seen  in  Figure 
4.9.  Note  also  the  magnitude  of  the  maneuver  in  this  solution:  Avz  ~  1.16  m/s. 


Right  Ascension 


Figure  4.6:  Maneuver  Fit  to  29  —  31  July,  tm  =  79477 
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Figure  4.7:  RA  vs  Dec  for  29  —  31  July,  tm  =  79477 


Right  Ascension 


Figure  4.8:  Maneuver  Fit  to  29  —  31  July,  tm  =  125340 
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Figure  4.9:  RA  vs  Dec  for  29  —  31  July,  tm  =  125340 

The  solution  for  the  final  convergence  value  tm  =  168380  is  shown  in  Figure  4.10. 
This  solution  is  unlikely  due  to  the  growing  separation  in  the  right  ascension  plot.  This 
growth  can  also  be  seen  in  the  right  ascension  versus  declination  plot  shown  in  Figure 
4.11.  Note  that  the  solution  does  not  give  a  good  fit  for  the  first  set  of  observations 
and  completely  misses  the  second  set  of  observations.  Note  also  the  large  magnitude 
of  the  maneuver  in  two  directions:  Ary  ~  0.72  and  Avz  ~  1.16  m/s.  In  many  cases, 
maneuver  magnitude  can  be  a  useful  discrimination  tool.  If  information  about  a 
particular  collocated  satellite  is  known  (such  as  thruster  type,  amount  of  fuel,  etc.) 
it  is  possible  to  rule  out  any  solutions  that  would  exceed  operational  capabilities. 

The  issue  of  different  solutions  for  different  initial  tm  values,  producing  apparent 
local  minima,  is  an  important  issue  to  investigate.  One  possible  explanation  would 
be  noise  inherent  in  the  data.  While  this  could  have  contributed,  the  fact  that  this 
same  phenomenon  occurred  in  the  simulated  real  data  (using  perfect  observations), 
presented  in  Section  4. 2. 2.1,  suggests  that  noise  is  not  the  driving  factor.  Another 
possible  explanation  could  be  that  more  than  one  maneuver  has  taken  place.  The 
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plausibility  of  this  situation  is  validated  in  another  data  set  as  presented  in  Section 
4.4.2.  Finally,  the  most  likely  explanation  lies  in  the  quantity  and  non-continuous 
nature  of  the  observations  available.  As  shown  by  the  maneuver  simulation  results, 
observability  could  prove  to  be  the  largest  obstacle  in  maneuver  estimation.  Given 
the  available  solutions  and  analysis  in  the  section  above,  however,  tm  =  79477  was 
determined  to  be  the  most  likely  solution. 


Right  Ascension 


Figure  4.10:  Maneuver  Fit  to  29  —  31  July,  tm  =  168380 


4-3.2  Data  Set  2:  30  July  -  1  August.  The  second  set  of  data,  consisting  of 
observations  obtained  on  30  July  -  1  August,  was  used  to  examine  a  possible  maneuver 
performed  by  Sat  2  between  the  31st  of  July  and  the  1st  of  August. 

As  seen  in  Figure  4.12,  the  solution  produced  by  the  linear  non-maneuver  model 
does  not  give  an  accurate  fit  to  these  data  arcs  either.  The  declination  prediction 
seems  to  follow  the  trend  in  the  data;  however,  the  right  ascension  does  not  fit  the 
data  arcs  for  31  July  or  1  August.  The  state  vector  solution  Xo  is  shown  in  Table 
4.10. 
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Figure  4.11:  RA  vs  Dec  for  29  —  31  July,  tm  =  168380 
Table  4.10:  State  Vector  Solution  for  Non- maneuver  Fit  to  30  July  -  1  August 


X0 

Solution 

dr 

1478.5 

r0S6 

-45779 

dz 

1830.2 

dr 

-0.17368 

r066 

-0.19512 

dz 

-0.68939 

The  data  arcs  were  then  run  in  the  maneuver  model,  and  as  mentioned  above, 
the  least  squares  algorithm  was  looped  through  for  each  observation  time. 

While  the  model  produced  numerous  possible  solutions,  no  convergence  was 
achieved  with  this  data  set.  Referring  back  to  Table  4.1,  the  data  obtained  on  30 
July,  the  first  day  of  this  data  set,  consisted  of  only  nine  observations  in  two  arcs 
separated  by  12  minutes.  It  is  therefore  most  probable  that  there  was  not  enough 
data,  particularly  on  the  first  day,  for  the  estimator  to  converge  upon  a  solution. 
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Figure  4.12:  Non-maneuver  Fit  to  31  July  -  1  August 

4-4  Raven  Data  -  2004 

Data  for  the  same  collocation  was  obtained  for  3  -  5  and  9  June  of  2004.  Two 
maneuvers  were  suspected;  one  maneuver  was  suspected  sometime  between  3  and 
4  June  and  another  sometime  between  5  and  9  June.  Due  to  the  unobservability 
between  5  and  9  June,  it  was  decided  to  focus  primarily  on  the  first  three  days  of 
observations.  Similar  to  the  process  used  for  29  July  -  1  August  2003,  this  data  was 
broken  into  two  different  runs.  The  first  data  set  consisted  of  3  and  4  June,  with 
the  second  data  set  consisting  of  observations  from  4  and  5  June.  These  three  Hies 
have  approximately  four  times  the  number  of  observations  than  those  from  July  2003. 
Table  4.11  gives  the  approximate  time  (given  in  hour  and  minute  of  that  day)  and 
length  (shown  in  the  number  of  observations)  of  each  arc  for  3  June  as  well  as  the  time 
between  observations  within  the  arc,  At.  The  data  set  for  4  June  is  similar  to  that 
of  3  June.  It  consists  of  13  data  arcs  starting  at  hour  7  and  going  through  hour  11. 
There  are  10  arcs  having  5  observations  each  and  three  arcs  having  4  observations. 
Only  two  arcs  were  obtained  for  hour  9,  however,  one  at  the  beginning  of  the  hour 
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Table  4.11:  Data  Arc  Time  and  Length  for  3  June 


June 

Hour 

Minute 

Obs 

At  (s) 

3 

9 

41-42 

5 

16 

10 

2-3 

5 

16 

10 

34-35 

5 

16 

10 

45-46 

5 

16 

10 

56-57 

5 

16 

11 

6-7 

5 

16 

11 

17-18 

5 

16 

11 

28-29 

5 

16 

11 

38-39 

4 

16 

11 

49-50 

5 

16 

12 

0-1 

5 

16 

12 

10-11 

5 

16 

12 

21-22 

5 

16 

12 

32-33 

5 

16 

12 

42-43 

5 

16 

12 

53-54 

3 

16 

13 

4-5 

5 

16 

13 

14-16 

5 

16 

13 

25-26 

3 

16 

and  one  at  the  end  of  the  hour.  This  break  in  data  can  be  observed  in  the  plots  to 
follow.  The  data  set  for  5  June  consists  of  18  arcs  having  5  observations  and  one  arc 
with  4  observations  starting  at  hour  9  and  going  through  hour  11. 

4-4-1  Data  Set  1:  3-4  June.  Using  the  same  procedure  as  that  described 
in  Section  4.3,  the  first  step  was  to  run  these  data  arcs  through  the  non-maneuver 
model.  The  plots  showing  right  ascension  and  declination  versus  time  look  as  though 
the  model  produced  an  accurate  solution.  See  Figure  4.13.  The  right  ascension  versus 
declination  plot,  however,  reveals  a  slight  discrepancy  in  declination,  particularly  in 
the  last  set  of  arcs.  See  Figure  4.14. 

The  data  was  then  run  through  the  maneuver  model  using  the  same  method  as 
stated  above.  Only  one  maneuver  time  and  state  vector  completely  converged.  An 
initial  maneuver  time  guess  between  approximately  7  hours  and  38  min  before  and  4 
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hours  after  t  =  80008  s  resulted  in  convergence  on  the  state  vector  shown  in  Table 
4.12  with  a  maneuver  time  of  tm  =  80008  (hour  7,  minute  55  on  4  June).  Note  that 

Table  4.12:  State  Vector  Solution  for  Maneuver  Fit  to  3-4  June 


X0 

tm  =  80008 

Sr 

7229.6 

r06d 

-1.0327e+005 

Sz 

-25696 

Sr 

0.49347 

r0Sd 

-1.095 

Si 

-2.1298 

Avr 

-0.03253 

Avg 

-0.0071178 

<1 

0.20463 

VjTl 

4.6267e-010 

the  main  component  of  the  maneuver  in  the  state  vector  solution  is  in  the  Sz,  or  n-s 
direction.  This  is  consistent  with  a  correction  in  declination. 

Figure  4.15  shows  right  ascension  and  declination  versus  time.  A  slightly  better 
fit  to  the  declination  can  be  seen  in  the  lower  plot.  A  significant  improvement,  how¬ 
ever,  is  discernable  in  the  right  ascension  versus  declination  plot.  Figure  4.16  confirms 
the  primary  burn  in  the  n-s  direction  results  in  a  better  fit  to  the  second  set  of  data 
arcs  from  4  June.  While  tm  =  80008  was  the  only  value  converged  upon  for  multiple 
different  values  of  initial  maneuver  time  guesses,  it  should  be  noted  that  numerous 
single  state  vector  and  maneuver  time  solutions  produced  reasonable  answers.  With¬ 
out  convergence,  however,  there  is  little  confidence  that  these  answers  are  legitimate 
solutions. 
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Right  Ascension 


Figure  4.13:  Non-maneuver  Fit  to  3-4  June 


Right  Ascension  Difference  (micro  rad) 


Figure  4.14:  Non-maneuver  RA  vs  Dec  for  3-4  June 
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Figure  4.15:  Maneuver  Fit  to  3  —  4  June,  tm  =  80008 


Right  Ascension  Difference  (micro  rad) 


Figure  4.16:  RA  vs  Dec  for  3  —  4  June,  tm  =  80008 
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4-4-2  Data  Set  2:  4  ~  5  June.  Because  the  maneuver  was  suspected  between 
3  and  4  June,  not  4  and  5  June,  it  would  be  expected  that  the  solution  from  the 
previous  section  would  fit  the  observations  from  5  June  if  propagated  forward  for  an 
additional  day.  As  shown  in  Figure  4.17,  this  assumption  was  not  validated.  The 
solution  for  3  and  4  June  did  not  fit  the  observations  for  the  5th.  Note,  however,  that 


Figure  4.17:  RA  vs  Dec  Solution  for  3  —  4  June  Propagated  to  5  June 

the  change  in  declination  from  4  to  5  June  (shown  by  an  upward  shift  from  the  second 
ellipse  to  the  new  observations)  is  equal  to  the  change  in  declination  from  3  to  4  June 
(shown  by  an  upward  shift  from  the  first  ellipse  to  the  second  ellipse).  This  would 
seem  to  indicate  another  maneuver. 

Due  to  this  discrepancy,  it  was  necessary  to  revert  back  to  running  the  data 
from  4-5  June  in  the  maneuver  model.  Similar  to  the  previous  results,  one  maneuver 
time,  tm  =  87311,  and  state  vector  solution  was  converged  upon.  This  solution  can 
be  seen  in  Table  4.13. 

As  seen  in  Figure  4.18,  not  much  useful  information  can  be  extracted  from  the 
plot  showing  right  ascension  and  declination  versus  time.  Right  ascension  plotted 
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Table  4.13: 


State  Vector  Solution  for  Maneuver  Fit  to  4-5  June 


X0 

tm  =  87311 

Sr 

3703.9 

r056 

-87686 

Sz 

-10936 

Sr 

0.66135 

r0S6 

-0.57611 

Si 

-2.5168 

Avr 

-0.015592 

A  vg 

0.005491 

< 

0.21955 

5tm 

1.0763e-011 

versus  declination,  however,  reveals  the  maneuver  model  solution  produced  a  more 
accurate  fit.  This  can  be  seen  in  Figure  4.19. 


Right  Ascension 


Figure  4.18:  Maneuver  Fit  to  4  —  5  June,  tm  =  87311 
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Figure  4.19:  RA  vs  Dec  for  4  —  5  June,  tm  =  87311 
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V.  Conclusions 


5. 1  Summary 

This  thesis  demonstrates  the  ability  to  estimate  maneuvers  of  a  collocated  satel¬ 
lite  in  geosynchronous  orbit  using  optical  observations  and  relative  orbit  determina¬ 
tion.  Based  on  the  research  conducted,  it  has  been  found  that  maneuver  estimation 
has  great  potential  to  become  a  valuable  tool  in  satellite  tracking  and  identification. 

5.2  Conclusions 

While  this  research  has  demonstrated  the  potential  of  maneuver  estimation,  it 
did  not  result  in  an  operationally  competent  model.  Solutions  were  obtained  for  three 
out  of  four  real  data  sets;  however,  several  issues  may  be  contributing  to  less  than 
confident  results. 

5.2.1  Observability.  Simulation  studies  with  perfect  data  show  that  the 
position,  velocity,  change  in  velocity,  and  time  of  maneuver  can  be  correctly  estimated 
with  a  set  of  observations  every  ten  minutes  over  the  course  of  an  orbit.  On  the  other 
hand,  current  amounts  of  data  obtained  for  tracking  satellites  have  been  found  to 
be  inadequate.  The  number  of  observations  obtained  for  the  first  data  set  (29  July 
-  1  August  2003)  generated  inconclusive  results.  The  second  data  set  (3-5  June 
2004),  however,  was  the  result  of  a  focused  effort  specifically  designed  to  obtain  more 
observations.  With  approximately  four  times  the  number  of  observations,  this  data 
set  yielded  promising  results. 

This  study,  therefore,  has  shown  that  if  maneuver  estimation  is  to  be  successful, 
a  concentrated  effort  must  be  employed  and  the  cluster  in  question  must  be  given 
a  higher  priority  so  as  to  obtain  more  observations  than  that  necessary  to  track  a 
regular  satellite.  The  level  of  observability  required  to  accurately  estimate  maneuvers 
in  GEO  is  far  greater  than  the  current  methodology  employed  to  simply  track  satellites 
or  clusters  of  satellites.  Obtaining  one  or  two  arcs  of  data  each  hour  for  three  to  four 
hours  has  been  shown  to  be  inadequate. 


5-1 


5.2.2  Higher  Order  Error  Sources.  Another  contribution  to  the  inconclusive 
maneuver  results  produced  with  running  the  real  data  in  the  maneuver  model  is  higher 
order  error  sources.  If  not  monitored,  small  errors  can  build  up  and  significantly  alter  a 
satellite’s  orbit,  thus  requiring  frequent  updates  to  the  orbit  determination  prediction. 
Combined  with  a  relatively  small  number  of  observations  over  the  course  of  an  orbit, 
error  sources  such  as  coordinate  system  errors,  solar  radiation  pressure,  and  other 
unmodeled  dynamics  could  serve  to  inhibit  convergence  within  the  model. 

5.2.3  Sequential  Maneuvers.  As  mentioned  in  Section  4.4.2,  sequential  ma¬ 
neuvers  are  also  a  possibility.  Due  to  the  recent  developments  in  propulsion  technol¬ 
ogy  many  companies  are  moving  towards  the  use  of  electric,  or  ion,  thrusters.  Ion 
thrusters  produce  a  higher  specific  impulse,  or  Isp,  than  their  chemical  counterparts, 
allowing  for  a  reduction  in  propellant  mass  and  an  increase  in  satellite  lifetime  [2], 
This  higher  Isp  usually  equates  to  lower  thrust  which  may  lead  to  an  increase  in 
the  number  of  burns  required  for  stationkeeping.  The  increased  efficiency  of  these 
thrusters,  however,  make  maneuvers  on  successive  days  a  possibility. 

5.3  Future  Work 

Given  the  conclusions  above,  there  are  several  areas  of  additional  research  which 
would  benefit  the  maneuver  estimation  process. 

One  such  area  is  to  study  the  effects  of  small  errors  to  include  solar  radiation 
pressure  and  other  unmodeled  dynamics,  coordinate  system  errors,  and  higher  order 
stationkeeping  perturbations.  By  systematically  isolating  different  potential  error 
sources  it  may  be  possible  to  quantify  the  effects  of  these  errors. 

Another  area  of  research  involves  determining  the  amount  of  data  required  to 
obtain  a  confident  solution.  This  would  also  include  evaluating  observation  quantities 
and  available  optical  assets  versus  necessary  optical  assets  across  the  globe. 

Given  the  increased  use  of  ion  thrusters,  the  possibility  of  multiple  maneuvers 
within  a  data  set  should  be  investigated. 
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Finally,  a  priori  information  on  actual  stationkeeping  maneuvers,  as  well  as 
satellite  capabilities,  would  greatly  assist  in  the  model  validation  process.  This  would 
involve  establishing  working  relationships  with  satellite  owners  and  operators. 
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